library(uds)

load("post.uds.RData")

load("post.decile.RData")
z.decile.95 <- apply(post.decile$z, 2, quantile, probs=c(.025, .975))
z.uds.95 <- apply(post.uds$z, 2, quantile, probs=c(.025, .975))
print(sum(z.decile.95[1,] > z.uds.95[2,] | 
          z.decile.95[2,] < z.uds.95[1,]))

z.decile.iqr <- apply(post.decile$z, 2, quantile, probs=c(.25, .75))
z.uds.iqr <- apply(post.uds$z, 2, quantile, probs=c(.25, .75))
print(sum(z.decile.iqr[1,] > z.uds.iqr[2,] |
          z.decile.iqr[2,] < z.uds.iqr[1,]))

sigmasq.decile.95 <- apply(post.decile$sigmasq, 2, quantile, probs=c(.025, .975))
sigmasq.uds.95 <- apply(post.uds$sigmasq, 2, quantile, probs=c(.025, .975))
print(sum(sigmasq.decile.95[1,] > sigmasq.uds.95[2,] | 
          sigmasq.decile.95[2,] < sigmasq.uds.95[1,]))

sigmasq.decile.iqr <- apply(post.decile$sigmasq, 2, quantile, probs=c(.25, .75))
sigmasq.uds.iqr <- apply(post.uds$sigmasq, 2, quantile, probs=c(.25, .75))
print(sum(sigmasq.decile.iqr[1,] > sigmasq.uds.iqr[2,] |
          sigmasq.decile.iqr[2,] < sigmasq.uds.iqr[1,]))

load("ppd.decile.RData")
ppd.hist(post.decile, ppd)

rm(post.decile, ppd)
load("post.vintile.RData")

z.vintile.95 <- apply(post.vintile$z, 2, quantile, probs=c(.025, .975))
z.uds.95 <- apply(post.uds$z, 2, quantile, probs=c(.025, .975))
print(sum(z.vintile.95[1,] > z.uds.95[2,] |
          z.vintile.95[2,] < z.uds.95[1,]))

z.vintile.iqr <- apply(post.vintile$z, 2, quantile, probs=c(.25, .75))
z.uds.iqr <- apply(post.uds$z, 2, quantile, probs=c(.25, .75))
print(sum(z.vintile.iqr[1,] > z.uds.iqr[2,] | 
          z.vintile.iqr[2,] < z.uds.iqr[1,]))

sigmasq.vintile.95 <- apply(post.vintile$sigmasq, 2, quantile, probs=c(.025, .975))
sigmasq.uds.95 <- apply(post.uds$sigmasq, 2, quantile, probs=c(.025, .975))
print(sum(sigmasq.vintile.95[1,] > sigmasq.uds.95[2,] |
          sigmasq.vintile.95[2,] < sigmasq.uds.95[1,]))

sigmasq.vintile.iqr <- apply(post.vintile$sigmasq, 2, quantile, probs=c(.25, .75))
sigmasq.uds.iqr <- apply(post.uds$sigmasq, 2, quantile, probs=c(.25, .75))
print(sum(sigmasq.vintile.iqr[1,] > sigmasq.uds.iqr[2,] | 
          sigmasq.vintile.iqr[2,] < sigmasq.uds.iqr[1,]))

load("ppd.vintile.RData")
ppd.hist(post.vintile, ppd)

rm(list=ls())
